module efso_mod

  use flogger
  use fiona
  use container
  use atm_formula_mod
  use data_translators_utils_mod, only: divide
  use const_mod, only: ens_ibeg, ens_iend, Rd, cp
  use math_mod
  use efso_utils_mod
  use namelist_mod
  use obs_data_mod
  use latlon_data_mod
  use grapes_reader_mod
  use era5_reader_mod
  use era5_regrid_mod

  implicit none

  private

  public efso_init
  public efso_run
  public efso_final

  type(timedelta_type) half_window_size

  type(grapes_reader_type) model_reader
  type(latlon_data_type), allocatable :: model_data_t0(:)
  type(latlon_data_type), allocatable :: model_data_tn(:)

  type(latlon_data_type) model_mean_data_t0
  type(latlon_data_type) model_mean_data_tn
  type(latlon_data_type) ref_data_tn
  type(latlon_data_type) fcst_err_tn

  real(r8), allocatable :: obs_impact(:)

contains

  subroutine efso_init(namelist_file)

    character(*), intent(in) :: namelist_file

    integer iens

    call process_init()
    call fiona_init()
    call namelist_init(namelist_file)

    half_window_size = create_timedelta(hours=3)

    call obs_data_init()

    ! Initialize model data.
    call model_reader%init(model_files_t0(1))
    ! NOTE: If we choose to parallel with ensemble members, we do not decompose domain.
    call proc%init(parallel_domain=.not. parallel_ensemble)
    if (parallel_ensemble) then
      ! Distribute ensemble members to different processes.
      call round_robin(ensemble_size, proc%np, proc%id, ens_ibeg, ens_iend)
      ens_ibeg = ens_ibeg + 1
      ens_iend = ens_iend + 1
    else
      ens_ibeg = 1
      ens_iend = ensemble_size
    end if
    allocate(model_data_t0(ens_ibeg:ens_iend))
    allocate(model_data_tn(ens_ibeg:ens_iend))
    do iens = ens_ibeg, ens_iend
      call model_data_t0(iens)%init(eid=iens)
      call model_data_tn(iens)%init(eid=iens)
    end do
    call model_mean_data_t0%init(eid=0)
    call model_mean_data_tn%init(eid=0)
    call ref_data_tn%init(eid=0)
    call fcst_err_tn%init(eid=0)

    if (proc%is_root_proc() .and. use_synop) then
      call log_error('The surface variables of initial ensemble members of GRAPES are all the same, cannot evaluate SYNOP.')
    end if

  end subroutine efso_init

  subroutine efso_run()

    integer iens, i, j

    ! Read observation data.
    call obs_data_read(model_reader%get_time() - half_window_size, model_reader%get_time() + half_window_size)

    ! Connect observation sites and model grids.
    call obs_data_interp_prepare()

    ! Read initial model data.
    do iens = ens_ibeg, ens_iend
      if (parallel_ensemble) then
        call log_notice('Open ' // trim(model_files_t0(iens)) // '.')
      end if
      call model_reader%open(model_files_t0(iens))
      call model_reader%get_var('orog', model_data_t0(iens)%zs)
      call model_reader%get_var('sp ', model_data_t0(iens)%ps )
      call model_reader%get_var('2t ', model_data_t0(iens)%t2 )
      call model_reader%get_var('2r ', model_data_t0(iens)%rh2)
      call model_reader%get_var('10u', model_data_t0(iens)%u10)
      call model_reader%get_var('10v', model_data_t0(iens)%v10)
      call model_reader%get_var('u  ', model_data_t0(iens)%u  )
      call model_reader%get_var('v  ', model_data_t0(iens)%v  )
      call model_reader%get_var('t  ', model_data_t0(iens)%t  )
      call model_reader%get_var('q  ', model_data_t0(iens)%q  )
      call model_reader%get_var('gh ', model_data_t0(iens)%z  )

      model_data_t0(iens)%ps = model_data_t0(iens)%ps * 0.01_r8 ! Pa -> hPa
      model_data_t0(iens)%t2 = model_data_t0(iens)%t2 - T0      ! K  -> degC
      model_data_t0(iens)%t  = model_data_t0(iens)%t  - T0      ! K  -> degC
      do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
        do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
          model_data_t0(iens)%q2(i,j) = specific_humidity_from_relative_humidity(model_data_t0(iens)%ps(i,j), &
                                                                                 model_data_t0(iens)%t2(i,j), &
                                                                                 model_data_t0(iens)%rh2(i,j))
        end do
      end do

      model_data_t0(iens)%time = model_reader%get_time()

      ! Interpolate model data to observation sites. Y
      call obs_data_interp_model_data(iens, model_data_t0(iens))
    end do
    ! Read forecast model data.
    do iens = ens_ibeg, ens_iend
      if (parallel_ensemble) then
        call log_notice('Open ' // trim(model_files_tn(iens)) // '.')
      end if
      call model_reader%open(model_files_tn(iens))
      call model_reader%get_var('orog', model_data_tn(iens)%zs)
      call model_reader%get_var('sp ', model_data_tn(iens)%ps )
      call model_reader%get_var('2t ', model_data_tn(iens)%t2 )
      call model_reader%get_var('2r ', model_data_tn(iens)%rh2)
      call model_reader%get_var('10u', model_data_tn(iens)%u10)
      call model_reader%get_var('10v', model_data_tn(iens)%v10)
      call model_reader%get_var('u  ', model_data_tn(iens)%u  )
      call model_reader%get_var('v  ', model_data_tn(iens)%v  )
      call model_reader%get_var('t  ', model_data_tn(iens)%t  )
      call model_reader%get_var('q  ', model_data_tn(iens)%q  )
      call model_reader%get_var('gh ', model_data_tn(iens)%z  )

      model_data_tn(iens)%ps = model_data_tn(iens)%ps * 0.01_r8 ! Pa -> hPa
      model_data_tn(iens)%t2 = model_data_tn(iens)%t2 - T0      ! K  -> degC
      model_data_tn(iens)%t  = model_data_tn(iens)%t  - T0      ! K  -> degC
      do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
        do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
          model_data_tn(iens)%q2(i,j) = specific_humidity_from_relative_humidity(model_data_tn(iens)%ps(i,j), &
                                                                                 model_data_tn(iens)%t2(i,j), &
                                                                                 model_data_tn(iens)%rh2(i,j))
        end do
      end do

      model_data_tn(iens)%time = model_reader%get_time()
    end do

    ! Calculate the ensemble mean of model. mean(X)
    call log_notice('Calculate ensemble mean.', pid=proc%id)
    call model_mean_data_t0%zero()
    call model_mean_data_tn%zero()
    do iens = ens_ibeg, ens_iend
      call model_mean_data_t0%add(model_data_t0(iens))
      call model_mean_data_tn%add(model_data_tn(iens))
    end do
    call model_mean_data_t0%global_sum()
    call model_mean_data_tn%global_sum()
    call model_mean_data_t0%divide(ensemble_size * 1.0_r8)
    call model_mean_data_tn%divide(ensemble_size * 1.0_r8)

    ! Interpolate ensemble mean to observation sites. H mean(X)
    call obs_data_interp_model_data(-1, model_mean_data_t0)

    ! Calculate the ensemble mean of simulated records on observation. mean(Y)
    call obs_data_calc_mean_simulated_records()

    ! Turn Y into Y'.
    call obs_data_calc_perb_simulated_records()

    ! Calculate innovation at each observation sites. mean(Y) - Yo
    call obs_data_calc_innov()

    ! Calculate ensemble perturbation (overwrite the data). X'
    call log_notice('Calculate ensemble perturbation.', pid=proc%id)
    do iens = ens_ibeg, ens_iend
      call model_data_t0(iens)%sub(model_mean_data_t0)
      call model_data_tn(iens)%sub(model_mean_data_tn)
    end do

    if (grapes_file_tn /= '') then
      call log_notice('Use ' // trim(grapes_file_tn) // ' as reference.', pid=proc%id)

      ! Read GRAPES analysis at tn as reference solution.
      call model_reader%open(grapes_file_tn)
      call model_reader%get_var('sp ', ref_data_tn%ps )
      call model_reader%get_var('2t ', ref_data_tn%t2 )
      call model_reader%get_var('2r ', ref_data_tn%rh2)
      call model_reader%get_var('10u', ref_data_tn%u10)
      call model_reader%get_var('10v', ref_data_tn%v10)
      call model_reader%get_var('u  ', ref_data_tn%u  )
      call model_reader%get_var('v  ', ref_data_tn%v  )
      call model_reader%get_var('t  ', ref_data_tn%t  )
      call model_reader%get_var('q  ', ref_data_tn%q  )
      call model_reader%get_var('gh ', ref_data_tn%z  )
    else if (era5_file_tn /= '') then
      call log_notice('Use ' // trim(era5_file_tn) // ' as reference.', pid=proc%id)

      ! Read ERA5 as reference solution.
      call era5_reader_run(era5_file_tn)

      ! Interpolate ERA5 onto model grids. Xt
      ref_data_tn%zs = model_data_t0(ens_ibeg)%zs
      call era5_regrid_run(ref_data_tn)
    end if

    ! Calculate forecast error.
    fcst_err_tn = model_mean_data_tn
    call fcst_err_tn%sub(ref_data_tn)

    call efso_core()

    call efso_write()

  end subroutine efso_run

  subroutine efso_final()

    call barrier()

    if (allocated(model_data_t0)) deallocate(model_data_t0)
    if (allocated(model_data_tn)) deallocate(model_data_tn)
    if (allocated(obs_impact   )) deallocate(obs_impact   )

    call era5_reader_final()
    call obs_data_final()
    call process_final()

  end subroutine efso_final

  subroutine efso_core()

    integer nens, nobs, ngrd_2d, ngrd_3d, ngrd
    integer iens, iobs, igrd, irec, imask, i, j, k

    real(r8), allocatable :: local_yb (:,:)
    real(r8), allocatable :: local_ybr(:,:)
    real(r8), allocatable :: yb (:,:)
    real(r8), allocatable :: ybr(:,:)
    real(r8), allocatable :: pa(:,:)
    real(r8), allocatable :: k_gain(:,:)
    real(r8), allocatable :: dy(:,:)
    real(r8), allocatable :: kdy(:,:)
    real(r8), allocatable :: local_xb_tn(:,:)
    real(r8), allocatable :: xb_tn(:,:)
    real(r8), allocatable :: xkdy(:,:)
    real(r8), allocatable :: xcxkdy(:,:)
    real(r8), allocatable :: djdy(:,:)

    real(r8), allocatable :: pa_inv(:,:)
    real(r8), allocatable :: eigvals(:)
    real(r8), allocatable :: eigvecs(:,:)

    real(r8), parameter :: tr = 280   ! K
    real(r8), parameter :: pr = 1.0e5 ! Pa
    real(r8), parameter :: wq = 1.0
    real(r8), parameter :: C_t = cp / tr
    real(r8), parameter :: C_q = lh**2 * wq / cp / tr
    real(r8), parameter :: C_p = Rd * tr / pr**2

    nens = ensemble_size
    ngrd_2d = proc%global_mesh%num_lon * proc%global_mesh%num_lat * 5
    ngrd_3d = proc%global_mesh%num_lon * proc%global_mesh%num_lat * proc%global_mesh%num_lev * 4
    ngrd = ngrd_2d + ngrd_3d
    nobs = (5 * obs_synop   %obs_size - count(obs_synop   %mask == 0)) + & ! u, v, t, q, p
           (4 * obs_raob    %obs_size - count(obs_raob    %mask == 0)) + & ! u, v, t, q
           (4 * obs_amdar   %obs_size - count(obs_amdar   %mask == 0)) + & ! u, v, t, q
           (2 * obs_profiler%obs_size - count(obs_profiler%mask == 0))     ! u, v

    ! Collect perturbated simulated records. Yb'
    allocate(local_yb (nobs,ens_ibeg:ens_iend))
    allocate(local_ybr(nobs,ens_ibeg:ens_iend))
    do iens = ens_ibeg, ens_iend
      i = 1
      associate (r => obs_synop%simulated_records, e => obs_synop%obs_err)
        imask = 1
        do irec = 1, obs_synop%records%size
          if (obs_synop%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%ua; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%u**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%va; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%v**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%ta; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%t**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%sh; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%q**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%p ; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%p**2; i = i + 1
          end if
          imask = imask + 1
        end do
      end associate
      j = 1
      associate (r => obs_raob%simulated_records, e => obs_raob%obs_err)
        imask = 1
        do irec = 1, obs_raob%records%size
          ! Mandatory levels
          do k = 1, obs_raob%simulated_records(irec,iens)%man%num_level
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%man%ua(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%u**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%man%va(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%v**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%man%ta(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%t**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%man%sh(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%q**2; i = i + 1
            end if
            imask = imask + 1
            j = j + 1
          end do
          ! Significant temperature levels
          do k = 1, obs_raob%simulated_records(irec,iens)%sigt%num_level
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigt%ua(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%u**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigt%va(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%v**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigt%ta(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%t**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigt%sh(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%q**2; i = i + 1
            end if
            imask = imask + 1
            j = j + 1
          end do
          ! Significant wind levels
          do k = 1, obs_raob%simulated_records(irec,iens)%sigw%num_level
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigw%ua(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%u**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigw%va(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%v**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigw%ta(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%t**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%sigw%sh(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%q**2; i = i + 1
            end if
            imask = imask + 1
            j = j + 1
          end do
          ! Tropopause levels
          do k = 1, obs_raob%simulated_records(irec,iens)%trop%num_level
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%trop%ua(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%u**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%trop%va(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%v**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%trop%ta(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%t**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%trop%sh(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%q**2; i = i + 1
            end if
            imask = imask + 1
            j = j + 1
          end do
        end do
      end associate
      associate (r => obs_amdar%simulated_records, e => obs_amdar%obs_err)
        imask = 1
        do irec = 1, obs_amdar%records%size
          if (obs_amdar%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%ua; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%u**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%va; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%v**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%ta; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%t**2; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            local_yb(i,iens) = r(irec,iens)%sh; local_ybr(i,iens) = local_yb(i,iens) / e(irec)%q**2; i = i + 1
          end if
          imask = imask + 1
        end do
      end associate
      j = 1
      associate (r => obs_profiler%simulated_records, e => obs_profiler%obs_err)
        imask = 1
        do irec = 1, obs_profiler%records%size
          do k = 1, obs_profiler%simulated_records(irec,iens)%pro%num_level
            if (obs_profiler%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%pro%ua(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%u**2; i = i + 1
            end if
            imask = imask + 1
            if (obs_profiler%mask(imask) == 1) then
              local_yb(i,iens) = r(irec,iens)%pro%va(k); local_ybr(i,iens) = local_yb(i,iens) / e(j)%v**2; i = i + 1
            end if
            imask = imask + 1
            j = j + 1
          end do
        end do
      end associate
    end do

    if (proc%is_root_proc()) then
      allocate(yb(nobs,nens), ybr(nobs,nens))
    else
      allocate(yb(1,1), ybr(1,1))
    end if
    call gather(proc%comm, ens_ibeg, ens_iend, local_yb , yb )
    call gather(proc%comm, ens_ibeg, ens_iend, local_ybr, ybr)
    deallocate(local_yb, local_ybr)

    ! Collect model data. Xb
    allocate(local_xb_tn(ngrd,ens_ibeg:ens_iend))
    do iens = ens_ibeg, ens_iend
      igrd = 1
      do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
        do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
          local_xb_tn(igrd,iens) = model_data_tn(iens)%u10(i,j); igrd = igrd + 1
          local_xb_tn(igrd,iens) = model_data_tn(iens)%v10(i,j); igrd = igrd + 1
          local_xb_tn(igrd,iens) = model_data_tn(iens)%t2 (i,j); igrd = igrd + 1
          local_xb_tn(igrd,iens) = model_data_tn(iens)%q2 (i,j); igrd = igrd + 1
          local_xb_tn(igrd,iens) = model_data_tn(iens)%ps (i,j); igrd = igrd + 1
        end do
      end do
      do k = proc%local_mesh%lev_ibeg, proc%local_mesh%lev_iend
        do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
          do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
            local_xb_tn(igrd,iens) = model_data_tn(iens)%u(i,j,k); igrd = igrd + 1
            local_xb_tn(igrd,iens) = model_data_tn(iens)%v(i,j,k); igrd = igrd + 1
            local_xb_tn(igrd,iens) = model_data_tn(iens)%t(i,j,k); igrd = igrd + 1
            local_xb_tn(igrd,iens) = model_data_tn(iens)%q(i,j,k); igrd = igrd + 1
          end do
        end do
      end do
    end do
    if (proc%is_root_proc()) then
      allocate(xb_tn(ngrd,nens))
    else
      allocate(xb_tn(1,1))
    end if
    call gather(proc%comm, ens_ibeg, ens_iend, local_xb_tn, xb_tn)
    deallocate(local_xb_tn)

    if (proc%is_root_proc()) then
      allocate(pa_inv(nens,nens))
      call matrix_mutiply('T', 'N', ybr, yb, pa_inv)
      do i = 1, nens
        pa_inv(i,i) = pa_inv(i,i) + nens - 1
      end do

      allocate(eigvals(nens))
      allocate(eigvecs(nens,nens))
      call log_notice('Calculate eigen values and vectors of matrix Pa.')
      !print *, sum(pa_inv)
      call symmetric_eigen_solver(pa_inv, eigvals, eigvecs)
      write(*, '(A20, A20)') 'Ensemble member', 'Eigen value'
      do j = 1, nens
        write(*, '(I20, F30.5)') j, eigvals(j)
        do i = 1, nens
          pa_inv(i,j) = eigvecs(i,j) / eigvals(j)
        end do
      end do
      allocate(pa(nens,nens))
      call matrix_mutiply('N', 'T', pa_inv, eigvecs, pa)
      deallocate(pa_inv, eigvals, eigvecs)

      allocate(k_gain(nens,nobs))
      call matrix_mutiply('N', 'T', pa, ybr, k_gain)

      ! Collect innovation. Yo - mean(Yb)
      allocate(dy(nobs,1))
      i = 1
      associate (innov => obs_synop%innov_records)
        imask = 1
        do irec = 1, obs_synop%obs_size
          if (obs_synop%mask(imask) == 1) then
            dy(i,1) = innov(irec)%ua; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            dy(i,1) = innov(irec)%va; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            dy(i,1) = innov(irec)%ta; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            dy(i,1) = innov(irec)%sh; i = i + 1
          end if
          imask = imask + 1
          if (obs_synop%mask(imask) == 1) then
            dy(i,1) = innov(irec)%p ; i = i + 1
          end if
          imask = imask + 1
        end do
      end associate
      associate (innov => obs_raob%innov_records)
        imask = 1
        do irec = 1, obs_raob%records%size
          ! Mandatory levels
          do k = 1, innov(irec)%man%num_level
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%man%ua(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%man%va(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%man%ta(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%man%sh(k); i = i + 1
            end if
            imask = imask + 1
          end do
          ! Significant temperature levels
          do k = 1, innov(irec)%sigt%num_level
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigt%ua(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigt%va(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigt%ta(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigt%sh(k); i = i + 1
            end if
            imask = imask + 1
          end do
          ! Significant wind levels
          do k = 1, innov(irec)%sigw%num_level
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigw%ua(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigw%va(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigw%ta(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%sigw%sh(k); i = i + 1
            end if
            imask = imask + 1
          end do
          ! Tropopause levels
          do k = 1, innov(irec)%trop%num_level
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%trop%ua(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%trop%va(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%trop%ta(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_raob%mask(imask) == 1) then
              dy(i,1) = innov(irec)%trop%sh(k); i = i + 1
            end if
            imask = imask + 1
          end do
        end do
      end associate
      associate (innov => obs_amdar%innov_records)
        imask = 1
        do irec = 1, obs_amdar%obs_size
          if (obs_amdar%mask(imask) == 1) then
            dy(i,1) = innov(irec)%ua; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            dy(i,1) = innov(irec)%va; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            dy(i,1) = innov(irec)%ta; i = i + 1
          end if
          imask = imask + 1
          if (obs_amdar%mask(imask) == 1) then
            dy(i,1) = innov(irec)%sh; i = i + 1
          end if
          imask = imask + 1
        end do
      end associate
      associate (innov => obs_profiler%innov_records)
        imask = 1
        do irec = 1, obs_profiler%records%size
          do k = 1, innov(irec)%pro%num_level
            if (obs_profiler%mask(imask) == 1) then
              dy(i,1) = innov(irec)%pro%ua(k); i = i + 1
            end if
            imask = imask + 1
            if (obs_profiler%mask(imask) == 1) then
              dy(i,1) = innov(irec)%pro%va(k); i = i + 1
            end if
            imask = imask + 1
          end do
        end do
      end associate

      allocate(kdy(nens,1))
      call matrix_mutiply('N', 'N', k_gain, dy, kdy)

      allocate(xkdy(ngrd,1))
      call matrix_mutiply('N', 'N', xb_tn, kdy, xkdy)

      igrd = 1
      do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
        do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
          xkdy(igrd,1) = xkdy(igrd,1) + 2 * fcst_err_tn%u10(i,j); igrd = igrd + 1
          xkdy(igrd,1) = xkdy(igrd,1) + 2 * fcst_err_tn%v10(i,j); igrd = igrd + 1
          xkdy(igrd,1) = C_t * (xkdy(igrd,1) + 2 * (fcst_err_tn%t2(i,j) + T0)); igrd = igrd + 1
          xkdy(igrd,1) = C_q * (xkdy(igrd,1) + 2 * fcst_err_tn%q2(i,j)); igrd = igrd + 1
          xkdy(igrd,1) = C_p * (xkdy(igrd,1) + 2 * (fcst_err_tn%ps(i,j) * 100)); igrd = igrd + 1
        end do
      end do
      do k = proc%local_mesh%lev_ibeg, proc%local_mesh%lev_iend
        do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
          do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
            xkdy(igrd,1) = xkdy(igrd,1) + 2 * fcst_err_tn%u(i,j,k); igrd = igrd + 1
            xkdy(igrd,1) = xkdy(igrd,1) + 2 * fcst_err_tn%v(i,j,k); igrd = igrd + 1
            xkdy(igrd,1) = C_t * (xkdy(igrd,1) + 2 * (fcst_err_tn%t(i,j,k) + T0)); igrd = igrd + 1
            xkdy(igrd,1) = C_q * (xkdy(igrd,1) + 2 * fcst_err_tn%q(i,j,k)); igrd = igrd + 1
          end do
        end do
      end do

      allocate(xcxkdy(nens,1))
      call matrix_mutiply('T', 'N', xb_tn, xkdy, xcxkdy)

      allocate(djdy(nobs,1))
      call matrix_mutiply('T', 'N', k_gain, xcxkdy, djdy)
      djdy = djdy

      allocate(obs_impact(nobs))
      do iobs = 1, nobs
        obs_impact(iobs) = 0.5_r8 * dy(iobs,1) * djdy(iobs,1)
      end do

      deallocate(pa, k_gain, dy, kdy, xkdy, xcxkdy, djdy)
    end if
    deallocate(yb, ybr, xb_tn)

    call barrier()

  end subroutine efso_core

  subroutine efso_write()

    type(linked_list_iterator_type) it
    integer imask, irec, i, k

    if (.not. allocated(obs_impact)) return
    if (.not. proc%is_root_proc()) return

    open(10, file=obs_impact_file)
    i = 1
    imask = 1; irec = 1
    it = linked_list_iterator(obs_synop%records)
    do while (.not. it%ended())
      select type (record => it%value)
      type is (synop_record_type)
        if (obs_synop%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'SYNOP', 'U', trim(record%station%name), &
            record%station%lon, record%station%lat, record%station%z, &
            obs_impact(i), record%ua, obs_synop%innov_records(irec)%ua
          i = i + 1
        end if
        imask = imask + 1
        if (obs_synop%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'SYNOP', 'V', trim(record%station%name), &
            record%station%lon, record%station%lat, record%station%z, &
            obs_impact(i), record%va, obs_synop%innov_records(irec)%va
          i = i + 1
        end if
        imask = imask + 1
        if (obs_synop%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'SYNOP', 'T', trim(record%station%name), &
            record%station%lon, record%station%lat, record%station%z, &
            obs_impact(i), record%ta, obs_synop%innov_records(irec)%ta
          i = i + 1
        end if
        imask = imask + 1
        if (obs_synop%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'SYNOP', 'Q', trim(record%station%name), &
            record%station%lon, record%station%lat, record%station%z, &
            obs_impact(i), record%sh, obs_synop%innov_records(irec)%sh
          i = i + 1
        end if
        imask = imask + 1
        if (obs_synop%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'SYNOP', 'P', trim(record%station%name), &
            record%station%lon, record%station%lat, record%station%z, &
            obs_impact(i), record%p , obs_synop%innov_records(irec)%p
          i = i + 1
        end if
        imask = imask + 1
        irec = irec + 1
      end select
      call it%next()
    end do
    imask = 1; irec = 1
    it = linked_list_iterator(obs_raob%records)
    do while (.not. it%ended())
      select type (record => it%value)
      type is (raob_record_type)
        ! Mandatory levels
        do k = 1, record%man%num_level
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'U', trim(record%station%name), &
              record%station%lon, record%station%lat, record%man%p(k), &
              obs_impact(i), record%man%ua(k), obs_raob%innov_records(irec)%man%ua(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'V', trim(record%station%name), &
              record%station%lon, record%station%lat, record%man%p(k), &
              obs_impact(i), record%man%va(k), obs_raob%innov_records(irec)%man%va(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'T', trim(record%station%name), &
              record%station%lon, record%station%lat, record%man%p(k), &
              obs_impact(i), record%man%ta(k), obs_raob%innov_records(irec)%man%ta(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'Q', trim(record%station%name), &
              record%station%lon, record%station%lat, record%man%p(k), &
              obs_impact(i), record%man%sh(k), obs_raob%innov_records(irec)%man%sh(k)
            i = i + 1
          end if
          imask = imask + 1
        end do
        ! Significant temperature levels
        do k = 1, record%sigt%num_level
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'U', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigt%p(k), &
              obs_impact(i), record%sigt%ua(k), obs_raob%innov_records(irec)%sigt%ua(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'V', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigt%p(k), &
              obs_impact(i), record%sigt%va(k), obs_raob%innov_records(irec)%sigt%va(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'T', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigt%p(k), &
              obs_impact(i), record%sigt%ta(k), obs_raob%innov_records(irec)%sigt%ta(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'Q', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigt%p(k), &
              obs_impact(i), record%sigt%sh(k), obs_raob%innov_records(irec)%sigt%sh(k)
            i = i + 1
          end if
          imask = imask + 1
        end do
        ! Significant wind levels
        do k = 1, record%sigw%num_level
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'U', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigw%p(k), &
              obs_impact(i), record%sigw%ua(k), obs_raob%innov_records(irec)%sigw%ua(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'V', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigw%p(k), &
              obs_impact(i), record%sigw%va(k), obs_raob%innov_records(irec)%sigw%va(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'T', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigw%p(k), &
              obs_impact(i), record%sigw%ta(k), obs_raob%innov_records(irec)%sigw%ta(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'Q', trim(record%station%name), &
              record%station%lon, record%station%lat, record%sigw%p(k), &
              obs_impact(i), record%sigw%sh(k), obs_raob%innov_records(irec)%sigw%sh(k)
            i = i + 1
          end if
          imask = imask + 1
        end do
        ! Tropopause levels
        do k = 1, record%trop%num_level
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'U', trim(record%station%name), &
              record%station%lon, record%station%lat, record%trop%p(k), &
              obs_impact(i), record%trop%ua(k), obs_raob%innov_records(irec)%trop%ua(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'V', trim(record%station%name), &
              record%station%lon, record%station%lat, record%trop%p(k), &
              obs_impact(i), record%trop%va(k), obs_raob%innov_records(irec)%trop%va(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'T', trim(record%station%name), &
              record%station%lon, record%station%lat, record%trop%p(k), &
              obs_impact(i), record%trop%ta(k), obs_raob%innov_records(irec)%trop%ta(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_raob%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'RAOB', 'Q', trim(record%station%name), &
              record%station%lon, record%station%lat, record%trop%p(k), &
              obs_impact(i), record%trop%sh(k), obs_raob%innov_records(irec)%trop%sh(k)
            i = i + 1
          end if
          imask = imask + 1
        end do
      end select
      irec = irec + 1
      call it%next()
    end do
    imask = 1; irec = 1
    it = linked_list_iterator(obs_amdar%records)
    do while (.not. it%ended())
      select type (record => it%value)
      type is (amdar_record_type)
        if (obs_amdar%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'AMDAR', 'U', trim(record%flight%name), record%lon, record%lat, record%p, &
            obs_impact(i), record%ua, obs_amdar%innov_records(irec)%ua
          i = i + 1
        end if
        imask = imask + 1
        if (obs_amdar%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'AMDAR', 'V', trim(record%flight%name), record%lon, record%lat, record%p, &
            obs_impact(i), record%va, obs_amdar%innov_records(irec)%va
          i = i + 1
        end if
        imask = imask + 1
        if (obs_amdar%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'AMDAR', 'T', trim(record%flight%name), record%lon, record%lat, record%p, &
            obs_impact(i), record%ta, obs_amdar%innov_records(irec)%ta
          i = i + 1
        end if
        imask = imask + 1
        if (obs_amdar%mask(imask) == 1) then
          write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'AMDAR', 'Q', trim(record%flight%name), record%lon, record%lat, record%p, &
            obs_impact(i), record%sh, obs_amdar%innov_records(irec)%sh
          i = i + 1
        end if
        imask = imask + 1
        irec = irec + 1
      end select
      call it%next()
    end do
    imask = 1; irec = 1
    it = linked_list_iterator(obs_profiler%records)
    do while (.not. it%ended())
      select type (record => it%value)
      type is (profiler_record_type)
        do k = 1, record%pro%num_level
          if (obs_profiler%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'PROFILER', 'U', trim(record%station%name), &
              record%station%lon, record%station%lat, record%pro%p(k), &
              obs_impact(i), record%pro%ua(k), obs_profiler%innov_records(irec)%pro%ua(k)
            i = i + 1
          end if
          imask = imask + 1
          if (obs_profiler%mask(imask) == 1) then
            write(10, '(3A10, 3F15.2, F20.2, 2F15.2)') 'PROFILER', 'V', trim(record%station%name), &
              record%station%lon, record%station%lat, record%pro%p(k), &
              obs_impact(i), record%pro%va(k), obs_profiler%innov_records(irec)%pro%va(k)
            i = i + 1
          end if
          imask = imask + 1
        end do
      end select
      irec = irec + 1
      call it%next()
    end do
    close(10)

    call log_notice('Observation impact results are written into ' // trim(obs_impact_file) // '.')

  end subroutine efso_write

end module efso_mod
